from shapely.geometry import shape, mapping
import geojson
import math
from geopy.distance import geodesic


def DegreeToRad(degree):
    return (degree /180)*math.pi

def GeoDist(point1,point2):
    """
    计算输入的两点之间的距离
    """
    radLat1=DegreeToRad(point1[1])
    radLat2=DegreeToRad(point2[1])
    delta_lon=DegreeToRad(point2[0]-point1[0])
    top_1=math.cos(radLat2)*math.sin(delta_lon)
    top_2=math.cos(radLat1)*math.sin(radLat2)-math.sin(radLat1)*math.cos(radLat2)*math.cos(delta_lon)
    top=math.sqrt(top_1*top_1+top_2*top_2)
    bottom=math.sin(radLat1)*math.sin(radLat2)+math.cos(radLat1)*math.cos(radLat2)*math.cos(delta_lon)
    delta_sigma=math.atan2(top,bottom)
    distance=delta_sigma*6378137

    return round(distance,3)


def GeoDist2(point1,point2):
    """
    使用geopy计算输入的两点之间的距离
    """
    lon1=point1[0]
    lat1=point1[1]
    lon2=point2[0]
    lat2=point2[1]
    distance=geodesic((lat1,lon1),(lat2,lon2)).meters
    return round(distance,3)



def GetVerticalDist(point1,point2,pointx):
    """
    计算输入的一点与弦之间的距离     

    输入：
    point1,point2 输入弦的两个端点
    pointx 待计算点弦距离的输入点

    输出：点弦距离（浮点型）
    """
    a=math.fabs(GeoDist(point1,point2))

    #当弦的两端重合时，点弦距离变为点间距离
    if a==0:
        return math.fabs(GeoDist(point1,pointx))
    
    #海伦公式求点弦距
    b=math.fabs(GeoDist(point1,pointx))
    c=math.fabs(GeoDist(point2,pointx))
    p=(a+b+c)/2
    s=math.sqrt(math.fabs(p*(p-a)*(p-b)*(p-c)))

    vertical_dist=s*2/a

    return vertical_dist


def DPCompress(input_points,threshold,depth_limit=100):
    """
    使用道格拉斯-普克法对矢量数据进行压缩

    input_points： 列表形式的输入点集，output_points： 列表形式的输出点集，
    threshold： 简化阈值
    depth_limit：递归深度，防止在点数过多时出现程序无响应
    """
    if len(input_points)<=2 or depth_limit==0:
        return input_points
    
    start_index=0
    end_index=len(input_points)-1
    index=start_index+1
    max_vertical_dist=0 #路径中离弦最远的距离
    key_point_index=0 #路径中的划分点

    while index<=end_index-1:
        #查找首末节点之间的节点到连线的距离
        cur_vertical_dist=GetVerticalDist(input_points[start_index],input_points[end_index],input_points[index])
        if cur_vertical_dist>max_vertical_dist:
            max_vertical_dist=cur_vertical_dist
            key_point_index=index  #将距连线最远的点设置为划分点
        index=index+1

        #将递归部分移至循环之外
        #递归划分路径
    if max_vertical_dist>=threshold:
        compressed_left=DPCompress(
            input_points[start_index:key_point_index+1],threshold,depth_limit-1
            )
        compressed_right=DPCompress(
            input_points[key_point_index:end_index+1],threshold,depth_limit-1
            )
        
        return compressed_left[:-1]+compressed_right   #去除重复的划分点
    else:
        return [input_points[start_index],input_points[end_index]]



def GetMeanErr(input_points,output_points):
    """
    计算压缩前后点位置的平均误差

    返回值：平均误差
    """
    err=0
    start_index=0
    end_index=len(output_points)-1
    i=0
    j=0
    input_points_new=[]
    output_points_new=[]

    for ipoint in input_points:
        ipoint=ipoint+[i]
        i=i+1
        input_points_new.append(ipoint)
    for opoint in output_points:
        opoint=opoint+[j]
        j=j+1
        output_points_new.append(opoint)


    while(start_index<=end_index-1):   #遍历所有关键点
        while(start_index+1<len(output_points)):
            point1_id=int(output_points_new[start_index][2])
            point2_id=int(output_points_new[start_index+1][2])

            id=point1_id+1
            # 找出point2_id在input_points中的位置
            for point in input_points_new:
                if point[2]>point1_id and point[2]!=input_points_new[-1][2]:
                    if point[0]==output_points_new[point2_id][0] and \
                            point[1]==output_points_new[point2_id][1]:
                        point3_id=point[2]

            while(id<=point3_id-1):
                err+=GetVerticalDist(output_points[start_index],
                                     output_points[start_index+1],input_points[id])
                id+=1

            start_index+=1

    return err/len(input_points)



def GeoJSONCompress(input_geojson, output_geojson,threshold):
    """
    使用道格拉斯-普克法压缩geojson数据

    input_geojson: geojson文件路径
    threshold: 阈值

    返回值: geojson文件路径
    """
    # Parse the input GeoJSON
    #geojson.load与geojson.loads函数有区别，load要求向其中传递一个TextIOWrapper
    #load则要求传递一个字符串
    input_points_num=0
    output_points_num=0
    input_points_sum=[]
    output_points_sum=[]
    with open(input_geojson, 'r',encoding='utf-8') as f:
        geojson_data = geojson.load(f)

    # Extract coordinates from the GeoJSON
    features = geojson_data['features']
    output_features = []

    for feature in features:
        geometry = feature['geometry']
        coordinates = geometry['coordinates']
        output_coordinates = []
        
        #geosjon文件中节点坐标是多层嵌套的，通过循环遍历每一个嵌套
        for nested_coords in coordinates:
            output_nested_coords=[]
            for coords in nested_coords:
                input_points=coords
                output_points=DPCompress(input_points, threshold)
                output_nested_coords.append(output_points)

                input_points_num+=len(input_points)
                output_points_num+=len(output_points)
                input_points_sum.extend(input_points)
                output_points_sum.extend(output_points)
            output_coordinates.append(output_nested_coords)


        output_feature = {'type': 'Feature',
                              'properties': feature['properties'],
                              'geometry': {'type': geometry['type'], 'coordinates': output_coordinates}}
        

        output_features.append(output_feature)
        print('input points:',input_points_num)
        print('output points:',output_points_num)
        
        #计算平均误差的函数可能存在问题
        # print('Mean Error:',GetMeanErr(input_points_sum,output_points_sum))


    transformed_geojson = {'type': 'FeatureCollection', 'features': output_features}
    with open(output_geojson, 'w',encoding='utf-8') as f:
        geojson.dump(transformed_geojson, f)
    print('done')


